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We consider the orbits of particles with spin in the Schwarzschild spacetime. Using the 
Papapetrou-Dixon equations of motion for spinning particles, we solve for the orbits and focus 
on those that exhibit chaos using both Poincare maps and Lyapunov exponents. In particular, we 
develop a method for comparing the Lyapunov exponents of chaotic orbits. We find chaotic orbits 
for smaller spin values than previously thought and find chaotic orbits with astrophysically relevant 
spin values. 
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I. INTRODUCTION 



A significant amount of effort has gone into gravita- 
tional wave detection over the last few decades. Many 
facilities, including LIGO, VIRGO, and GEO are dedi- 
cated to their detection and analysis. In addition, consid- 
erable effort has been made in modeling possible sources 
for gravitational waves. These include compact binaries 
with masses on the order of a few solar masses. These de- 
tectors are particulaly well tuned to such sources. How- 
ever, the proposed space based gravitational wave detec- 
tor, LISA, is better tuned to detecting signals such as 
those coming from compact objects orbiting supermas- 
sive black holes. 

Such extreme mass ratio inspirals (EMRI) [l[ have 
been studied extensively in recent years. In particular, 
in the test particle limit, the dynamics of these systems, 
together with the inclusion of spin, have been studied 
in a number of works [2 6]. A subset of these studies 
has considered the question of chaotic orbits in both the 
Schwarzschild [4] and Kerr [5, 6] spacetimes. This work 
has shown that chaotic orbits are possible in these sys- 
tems and is a consequence of the spin orbit coupling. 
However, these same studies suggest that chaotic orbits 
exist only when the orbiting particle has an unphysically 
large amount of spin. Nonetheless, the parameter space 
in which to search for chaos in these systems is large 
and has not been fully explored. It has been shown for 
both the Schwarzschild 0) and Kerr Q spacetimes that 
chaotic orbits change the character of the energy spec- 
trum of their gravitational waveform. Because chaotic 
orbits might also lead to chaotic gravitational wave sig- 
nals, and such signals may be more difficult to detect [jj, 
the more we can say about these systems the better. 

We consider the orbits of spinning test particles in 
Schwarzschild. Our approach to studying the possibil- 
ity of chaos in this system is to use Poincare sections as 
an indicator of chaos similar to Suzuki and Maeda [J]. 
Further, we use and extend a somewhat more sophisti- 
cated method of calculating the Lyapunov exponent of 



these orbits as employed by Hartl HQ. As part of this, 
we present a method for extracting an improved predic- 
tion of the Lyapunov exponent which allows us to dis- 
tinguish more carefully small Lyapunov exponents from 
zero. This analysis reveals some new classes of chaotic 
orbits. These chaotic orbits include a previously unsus- 
pected class of orbits that have spins for the inspiraling 
member that may be obtainable by astrophysical sys- 
tems. 

The remainder of the paper is constructed as follows. 
We describe the equations of motion next together with 
our choice of supplementary condition. In Sections IIIII 
and IIV1 we describe our methods for determining the 
chaos of various orbits, in Section[V]we present the results 
of our integrations and conclude in Section I VII 



II. THE PAPAPETROU-DIXON EQUATIONS 

To model a spinning test particle we use the equations 
of motion of Papapetrou [l(| and Dixon [TTJ] . These equa- 
tions describe a spinning particle in the pole-dipole ap- 
proximation. That is to say, they describe the particle 
as a mass monopole and spin dipole. These equations 
generalize geodesic motion to a spinning particle and are 
written in terms of the momentum, P a , and antisymmet- 
ric spin tensor, S ab , of the particle. They can be written 
as 
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where V a is the particle's velocity and, again, S ab defines 
the spin of the particle. As can be seen from Eq. ([2]) the 
momentum, P a , is not simply a rescaling of the velocity 
vector. Indeed, the momentum is defined as 



P a = f iV a - V b V c V c S : 



ab 



(3) 



where we will take fj, as the mass of our spinning test 
particle. 

As shown by Semerak Q this definition of the momen- 
tum can be combined with Eq. ^ to solve for V a is terms 
of the momentum and spin of the particle 
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As they stand, the equations of motions are underde- 
termined. This is a well known issue with these equa- 
tions and several supplementary conditions have been 
suggested in order to address this problem [Icil - [l2l ]. In 
this work, we choose to use the supplementary condition 
S ab P a — 0. This condition, in effect, picks out a center of 
mass frame for the particle. (For more discussion of this 
and other possible supplementary conditions see [T^ . [l3| .) 
Additionally, this supplementary condition implies that 
the spin tensor, S ab , has at most three independent com- 
ponents. As a result, it is possible to reformulate the 
equations of motion in terms of a spin vector, S a , which 
we will do below. 

As our background spacetime, we will choose spheri- 
cally symmetric Schwarzschild spacetime. This provides 
a number of conserved quantities with which to calculate 
test particle orbits, even with the assumed spin. In par- 
ticular, it can be shown that for a Killing vector X a the 
quantity 



C = X a P a + -s bc v b x c 



(5) 



is a constant of the motion. As a result, we can define 
the following two constants of the motion 
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In addition to these constants of the motion, the to- 
tal spin of the particle is also conserved. This quantity, 
hereafter S, is defined as the positive root of 



q2 ^~ Q QG,b 



(8) 



In order to understand the physical constraints on S, 
recall that lengths are measured in terms of the mass, m, 
of the central object of the spacetime and the momentum 
of the orbiting particle is measured in terms of its mass, 
/i. We might then think of S as being a unitless number 
multiplying m/i. Said another way, Sm/i = I where / is 
the spin angular momentum of the particle. Note that 
the Papapetrou equations are valid in the test particle 
approximation and therefore only hold for /i <g; m. Be- 
cause of this, the physical spin of the test particle must 
be much smaller than one in these units. This can be 
seen perhaps most clearly by considering the spinning 
test particle to be an extremal Kerr black hole orbiting 
around a supermassive black hole. For such an extremal 
black hole, its angular momentum is I = /i 2 . This leads 
to a total spin of 
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Using this argument, Hartl |5j estimates physical spins 
as being between about 10 -4 and 10~ 6 in S. 

For numerical simplicity we follow Suzuki and 
Maeda [J] and Hartl [a, Hi and modify the equation of 



motion by working with the spin vector. This quantity 
can be defined by 



(10) 



where £ a bcd is the totally antisymmetric tensor density. 
On making the following convenient definitions 
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the equations of motion [24j become 
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With this substitution, the velocity and constants of the 
motion are now defined by 
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Initial Conditions 



In order to characterize each orbit in as simple a way 
as possible, we make convenient choices in our initial con- 
ditions, working at all times with Schwarzschild coordi- 
nates. As every bound orbit will have turning points at 
which the radial velocity will be zero, we choose to be- 
gin all orbits with P r = 0. Also, because the spacetime 
is spherically symmetric we can set Pg — and begin 
motion in the equatorial plane. 

Under these conditions, our normalization, P a P a = 
— 1, and supplementary conditions, S a P a — 0, reduce to 
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respectively. These relations allow us to express the total 
spin as 
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FIG. 1: Three Poincare sections of the phase space trajec- 
tories of nonchaotic particle orbits. Each has L = 3.6 and 
from the innermost orbit out E\ — 0.9499, E2 = 0.9511, and 
E3 = 0.9534. Notice that the trajectories are confined to the 
surface of the tori which intersect the section. 



FIG. 2: A close up view of the orbits from Fig. [T] Notice 
that this zoomed in view continues to show the trajectories 
are constrained to the tori. 
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allows us to parametrize the spin vector components in 
terms of the angles a and /3. These angles are analogous 
to the 4> and of spherical polar coordinates respectively 
and have their origin at the test particle's center of mass. 

With these definitions we can specify an orbit by five 
initial conditions. The determining quantities are: r, P^, 
S, a, and f3. These correspond to the coordinate distance 
of the test particle from the central mass, its momentum 
in the </> direction, and the magnitude and orientation of 
its spin. 



III. MEASURING CHAOS 

In order to gain confidence in deciding whether a par- 
ticular orbit is chaotic or not we use two tests for chaos. 
The first is to check for the breaking up of the KAM tori 
in a Poincare section of the phase space [HI, [l6| . Follow- 
ing Suzuki we choose the section defined by the r — P r 
plane in the phase space, where r is the coordinate dis- 
tance the test particle is from the center of the central 
mass and P r is the conjugate momentum. In the case 
where the particle has no spin, typical sections look like 
ovals as in Figure [TJ A closer look at the these plots 
(Figure [2] ) shows that the phase space trajectories are 
confined to the surface of a torus. 

When considering sections for spinning test particle or- 
bits we look for this clean oval to break up as in Figure 
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FIG. 3: Three Poincare sections of the phase space trajecto- 
ries of chaotic particle orbits. From the inside orbit out the 
orbits have spin values Si = 1.515 S2 ~ 1.581 and S3 = 1.627. 
Notice that the trajectories stay close to the surface of the tori 
similar to Figure[TJ but do not have the thinly denned surface. 



[3] A closer look at the bands (Figure 2] ) clearly shows 
phase space trajectories have left the torus surface. The 
orbits that produce sections like this are close in to the 
black hole and are in agreement with the sections pro- 
duced by Q. 

The second method we use to look for chaos in this 
system is to calculate Lyapunov exponents. If one defines 
the distance D between two phase space trajectories as 



D (t) = d a e 



xt 



(24) 



where do is the initial separation between the two tra- 
jectories, then A is defined to be the Lyapunov expo- 
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FIG. 4: A close up view of the orbits from Fig. [3] Notice 
that this zoomed in view continues to show the trajectories 
are not constrained to the tori. Compare to Fig. [5] 



ncnt. When A is greater than zero the system is said to 
be chaotic. In Hamiltonian systems, of which the Papa- 
petrou equations are an example, the Lyapunov exponent 
can never be less than zero. A negative exponent would 
indicate some kind of attractor in the phase space, but 
these do not appear in conservative systems. 

In practice, calculating A is a challenge. The Jacobian 
method used by Hartl has the advantage of following just 
one phase space trajectory. The basic idea is that we can 
track the growth of a vector in the tangent space of the 
trajectory to compute the Lyapunov exponent (for more 
detail see @ and @). The equation that dictates the 
evolution of this vector is 



dr 



Df •£ 



(25) 



where £ is a vector in the tangent space, and Df is the 
Jacobian matrix of the system. 

To understand this method, consider a high dimen- 
sional ellipsoid in the phase space defined by some set 
of initial conditions. As the system evolves away from 
these initial conditions, this ellipsoid will become warped. 
However, by virtue of Liouville's theorem, the phase 
space volume of this ellipsoid will be conserved. As a 
result, any axis of the ellipsoid that corresponds to a 
chaotic coordinate will stretch. This will necessitate that 
at least one other axis will contract. The system can have 
more than one chaotic coordinate, but for each chaotic 
coordinate there must be a nonchaotic coordinate whose 
value converges at the same rate as the chaotic coordi- 
nate's value diverges. In other words, for every Lyapunov 
exponent that corresponds to a chaotic axis there is an 
exponent with the same magnitude but opposite sign. 

The Jacobian method finds the largest Lyapunov ex- 
ponent by evolving a vector, which is defined by the ini- 
tial conditions of the system, in the tangent space of the 



phase space. As the system evolves, this vector lines up 
with the direction of greatest stretching. By considering 
the magnitude of this vector as a function of the proper 
time r, we can then define the largest Lyapunov exponent 
of the system by 



A = lim 



(26) 



where £o denotes the initial tangent vector. As our nu- 
merical integrations cannot continue indefinitely, we will, 
in our calculations, refer to the Lyapunov exponent as a 
function of r, 



A(r) 



(27) 



where we have taken |£o| = 1- 

For this last analysis we have denoted the magnitude 
of a vector £ by |£|. Recall that this vector lives in the 
tangent space to the phase space of our physical system. 
This leads to uncertainty about what norm to use when 
calculating a vector magnitude. Eckmann [l7| shows that 
when calculating Lyapunov exponents, different norms 
may lead to different values, but the sign of the expo- 
nent will not be affected. With this in mind, we use the 
Euclidean norm for simplicity when calculating vector 
magnitudes. 

This method can be extended to find the Lyapunov ex- 
ponent corresponding to each axis of the stretching ellip- 
soid. Hartl implements this extended method @ in some 
cases and shows that the Lyapunov exponents do come 
in opposite sign pairs for the spinning particle system. 
His results also indicate that the direction of greatest 
stretching is not in the direction of any one coordinate 
or conjugate momentum. Thus, when we find the largest 
exponent we do not expect it to correspond to a partic- 
ular coordinate or momentum. 



IV. METHOD FOR COMPARING LYAPUNOV 
EXPONENTS 

Because the Jacobian method requires the Lyapunov 
exponent to be defined in terms of a limit we can in prac- 
tice only approximate its value. In previous work, a par- 
ticular orbit was allowed to evolve for some set amount 
of time and with the corresponding value of A(r) taken 
as the approximate exponent. 

A problem that arises with this method is that dif- 
ferent orbits converge to their Lyapunov exponents at 
different rates. In particular, the Lyapunov exponent for 
the case S — approaches zero much more slowly than 
for any other similar orbit with small, nonzero spin. To 
address this issue, we have developed a different method 
which both reduces computation time and predicts the 
Lyapunov exponent. 

Consider the plots of A(r) as shown in Figure [5] and 
Figure [5] These are typical examples of how A(r) con- 
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FIG. 5: We see the jagged Lyapunov function converging to 
zero for an orbit with S — 0. We also plot the curve fit to 
the data. This orbit begins at r = 6m with P^ = 3.6. The 
predicted Lyapunov exponent for this orbit is —1.68 x 10~ 5 
with an RMS error of is 1.23 x 10" 5 . 



verges for a zero spin orbit and a nonzero spin orbit re- 
spectively. Notice that the functions are modeled well by 
the fit 

f(r) = ax + £ (28) 

where a\, a2, and 03 are constants that are varied until 
the root mean square error between the fit and A(r) is 
minimized. More explicitly, the constants are varied to 
minimize 

N , v 2 

err(T) = Y,(ai + ^-\(T i )) . (29) 

It is then easy to define the RMS error for the fit as 



RMS = 




The fitting function sets the Lyapunov exponent of the 
system to be a\. Because our system is conservative we 
cannot have negative Lyapunov exponents. However, in 
the orbit corresponding to Figure [5] the root mean square 
error of the fit is 1.23x 10^ 5 while the calculated exponent 
is nonzero and just a bit outside this error range. Based 
on other results, such as Figure||)l as well as our method's 
consistency with the results of [I| and [f|, we are led to 
belive that the fit slightly underestimates the Lyapunov 
exponent. Another example is the orbit corresponding 
to Figure[6]in which the RMS error is 6.1 x 10~ 5 keeping 
zero well outside the error bars of the fit. 

This model of the Lyapunov function also gives a mea- 
sure of how quickly the exponent converges. In work- 
ing with this model we found that orbits with large spin 
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FIG. 6: We see the jagged Lyapunov function converging to a 
nonzero value for an orbit with S — 0.5. The initial conditions 
are r = 6m, — 3.6, a = j, and /3 = We also plot the 
curve fit to the data. The predicted Lyapunov exponent is 
3.787 x 10~ 4 with well outside the RMS error, 6.1 x 10~ 5 , 
of the fit. 



had Lyapunov functions that, in general, converged faster 
than those with small spin. It is possible that in previ- 
ous work similar Lyapunov exponents with slower con- 
vergence rates might have been discounted as numerical 
error. The difference in convergence rates between the 
zero point and these small spin exponents is much smaller 
than the difference between the zero point and the high 
spin exponents. If these differing rates of convergence 
are not taken into account, terminating each orbit after 
some predetermined time seems natural. Unfortunately, 
our experience shows that this can innate the value of 
the Lyapunov exponent in the zero spin case. Because 
that case has the slowest rate of convergence, it can then 
be difficult to distinguish it from cases with nonzero spin 
and nonzero Lyapunov exponents. 

We avoid this problem by comparing Lyapunov expo- 
nents once a a uniform degree of convergence has been 
reached rather than a particular time. While evolving the 
system we fit the Lyaponv function to Eq. (|28[) and con- 
tinue the evolution until the derivative of the fitting func- 
tion reaches a predefined tolerance close to zero. When 
the derivative of the fit has become sufficiently small, 
specifically when the value is on the order of 10~ 8 , we 
say that the function has converged. In this way we 
compare exponents which have all converged the same 
amount rather than comparing by the overall time of 
evolution. We find that this method reveals more of the 
chaotic nature of these orbits. 

Further, we can use the fit curve to predict the Lya- 
punov exponent for a given orbit without integrating for 
infinite time. These predicted values, together with an 
estimate of their error, give us confidence that the true 
exponent is within the corresponding range of values and, 
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FIG. 7: These Lyapunov exponents are given by the constant 
term in the curve fitting model, Eq. (|28[) . with root mean 
square error bars. The initial spin orientation is a = ? and 
P = j and constants of the motion for the spinless case are 
L = 3.6 and E = 0.9522. Compare with Figures [8] and \W\ 
Note that the orbits transition from nonchaotic to chaotic 
around a spin of S = 0.5. 

more importantly, provide a clear differentiation between 
positive and zero values. 

V. RESULTS 

A primary motivation for our work is to determine 
whether or not astrophysically relevant spins can give 
rise to chaotic orbits. We first consider orbits similar 
to those considered by Suzuki and Maeda and then look 
at other types of orbits. Because the orbits considered 
by Suzuki and Maeda are confined to high curvature re- 
gions, the additional orbital types we consider are those 
that remain in regions of low curvature and orbits that 
traverse both high and low curvature regions. For each 
orbital type we begin with a set of initial conditions with 
5 = and calculate the Lyapunov exponent. We then in- 
crease the magnitude of the spin and plot the Lyapunov 
exponents for each spin. 

A. High Curvature Orbits 

We first consider orbits similar to those used by Suzuki 
and Maeda |4j . These orbits remain close to the center of 
the spacetime throughout the system's evolution. This 
confines the orbits to regions of high curvature which 
allows the coupling between the curvature and spin to 
have a large effect on the particle's motion. 

One difference between the orbits we consider and 
those used by Suzuki and Maeda is the orbital angular 
momentum of the spinning particle. They consider orbits 



FIG. 8: These Lyapunov exponents are given by the constant 
value in the curve fitting model with the root mean square 
value providing the error bars. The initial spin orientation is 
P — tt and constants of the motion for the spinless case are 
L = 3.6 and E = 0.9522. All spin values higher than S = 0.3 
cause the particle to cross the event horizon. Note that lower 
spin values give exponents of the same order of magnitude as 
those in Figure [7] 



with L = 4.0 whereas we begin our analysis with orbits 
having angular momentum of L = 3.6. Our initial values 
for energy E — 0.9522 and an initial radius of r = 6m are 
chosen to make the orbit close to circular. We also keep 
track of the initial conditions in spin that produce these 
orbits. In particular, we consider what spin magnitudes 
produce chaotic orbits and how the spin vector's initial 
orientation affects the Lyapunov exponents. In the first 
case the spin vector's initial orientation is given by a = ? 
and ft = f . 

As Figure [7] shows, the Lyapunov values for this orbit 
can be put into one of two groups. From the maximal 
spin value of one down to about S = 0.5 there seems 
to be a well defined trend with nonzero Lyapunov expo- 
nents and corresponding chaotic orbits. These values are 
in agreement with those reported by Suzuki and Maeda 
as well as Hartl. From S = 0.5 down to zero spin the or- 
bits are not chaotic. Again, we find excellent agreement 
with 0] and @. 

The initial spin vector orientation can affect the mag- 
nitude of the Lyapunov exponent. As an example, if 
we change the spin orientation but otherwise keep the 
same initial conditions as before, the particle's behavior 
changes. When /3 = n, the spin is pointed down perpen- 
dicular to the equatorial plane, the spin orbit coupling 
causes the particle to be pulled in closer to the center 
of the spacetime. We find that when the spin angular 
momentum and orbital angular momentum are parallel 
there is a repulsive spin orbit interaction and when the 
momenta are antiparallel there is an attractive interac- 
tion. These results agree with Wald's 18] analysis of 
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FIG. 9: These Lyapunov exponents correspond to the small 
spin case of Figure [8] with f3 — n. Note the monotonically 
decreasing, yet still positive, values of the exponent for de- 
creasing spin values. This trend continues for spins as low 
as S — 0.015. However, this trend would appear to be lost 
for the very smallest values of spin. Indeed, we feel confident 
to treat these values as zero with correspondingly nonchaotic 
orbits. 



similar systems. When this interaction becomes strong 
enough the particle will no longer exhibit a bound orbit, 
and is captured by the black hole. In Figure [5] we see 
the Lyapunov exponent values for increasing spin values. 
Notice that after S » 0.3 there is no data. These data 
points are not included because for S > 0.3 the particle 
is captured by the black hole. 

Notice also that the exponents in Figure [5] appear to be 
nonzero below S = 0.05. Since the interaction between 
the spin angular momentum and the orbital angular mo- 
mentum pulls the particle in closer, the particle traverses 
higher curvature regions of the spacetime than the orbits 
described by Fig [71 This spin orientation is the only 
difference between the two cases. Thus, on comparing 
Figure [7] and Figure [5] we can see that the orientation of 
the spin can have a dramatic effect on the dynamics of 
the particle. 

In Figure [5] we zoom in on the small spin value orbits. 
This data clearly shows positive Lyapunov exponents for 
spins as small as S ~ 0.015. This lower bound in spin is 
considerably less than the bounds given in [3| and @. 
Notice that exponents in Figure[5]are small in comparison 
to Figure [7J Using our particular fitting method here to 
predict the Lyapunov exponent was crucial in being ablt 
to distinguish these small exponents from zero. 

When we set the initial spin orientation to (3 — 0, as 
in Figure [TU1 the particle is pushed further out from the 
center instead of pulled closer in. In this case the spin 
vector and orbital angular momentum vector are paral- 
lel. In this case the particle is never captured by the 
black hole. Instead the particle's orbit stays in areas of 



FIG. 10: These Lyapunov exponents are given by the constant 
value in the curve fitting model with the root mean square 
value providing the error bars. The initial spin orientation is 
P = and constants of the motion for the spinless case are 
L — 3.6 and E = 0.9522. No large exponents (compare to 
Figure 0) appear in this configuration. 



the spacetime with slightly less curvature. Also, the Lya- 
punov exponents stay at the same order of magnitude as 
the smaller values from Figure [7] for all spin values, but 
unlike Figure [7] no nonzero spins correspond to a zero 
exponent. In Figure QT] we look more closely at the small 
spin values for the same orbit. We notice the mono- 
tinicity in the predicted exponent values, and a nonzero 
exponent for spin values as low as S = 0.01. 

We also consider orbits which remain in areas of low 
curvature. Choosing nearly circular orbits at radii of 
100m and calculating the Lyapunov exponents as above, 
we find such orbits to have a uniformly zero Lyapunov 
exponent for every spin orientation which we consider 
above. When the radii is sufficiently reduced, r w 19m, 
positive exponents are again obtained. These positive ex- 
ponents first occur for very large spin, but as the radius 
continues to decrease, less spin is required to achieve a 
chaotic orbit. 



B. Knife Edge Orbits 

Knife edge orbits are some of the more interesting 
test particle trajectories allowed in black hole space- 
times. These orbits have large scale precession and exe- 
cute small tight loops around the center of the system. 
Their dynamics can be understood by considering the 
effective potential in the Schwarzschild spacetime. For 
large enough angular momentum this potential has a 
sharp peak close to the center of the spacetime. These 
"knife edge" potentials are what give these orbits their 
particular dynamics. [25j (For more discussion of these 
orbits see |2(|.) These orbits are also known as "zoom 
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FIG. 11: These Lyapunov exponents are given by the constant 
value in the curve fitting model with the root mean square 
providing the error bars. The initial spin orientation is /? = 
and constants of the motion for the spinless case are L — 3.6 
and E = 0.9522. Note that unlike Figure [7| the exponents all 
appear to be nonzero. No large exponents (compare to Figure 
appear in this configuration. Note the continuous nature 
to the exponents and the apparently positive value for spin 
as low as S = 0.01. 



whirl" orbits; the name being used to describe their be- 
havior, namely they zoom in from large radius and whirl 
about the center. (For a comprehensive cataloging of 
these orbits see (2lj): 

As we consider knife edge orbits for spinning test par- 
ticles we focus on the geometry surrounding their path. 
Unlike the orbits we have considered so far, these move 
through both high and low curvature regions of the space- 
time. These orbits start out relatively far from the center 
of the spacetime, where the curvature is comparatively 
low. But during the course of their orbits they execute 
several small orbits in much higher curvature regions. 
These types of orbits help us determine whether chaotic 
orbits must remain in high curvature or just pass through 
them regularly. 

In Figure [12] we plot the Lyapunov exponents for a 
knife edge orbit that has an outer radius of r = 100m and 
makes three small radius loops for each large radius orbit. 
The orientation of the spin is /3 = 0. The angular mo- 
mentum of the particle is L — 3.9246 and which is nearly 
the angular momentum of the particles investigated by 
Suzuki and Maeda. We choose this spin orientation to 
keep the particle from being captured by the black hole. 
Recall that this orientation gives an effective centrifugal 
force which pushes the particle away from the black hole. 
This keeps the particle from being captured, but also re- 
duces the number of inner loops traversed in each orbital 
period. Thus, as the spin increases, the particle can be 
thought to be retreating from the knife's edge. 

One important aspect to notice is that the S = orbit 



FIG. 12: These Lyapunov exponents correspond to knife edge 
orbits. Notice that the S = case is very nonzero. This is 
referred to as a "chaos mimic." The spinless orbit begins at 
r = 100m and has three inner loops at small radius for each 
large scale orbit. The initial spin orientation of ft = 0. In the 
spinless case E = 0.997 and L = 3.9246. 
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FIG. 13: These Lyapunov exponents correspond to the same 
knife edge orbits as Figure [121 Notice that zero is outside the 
error bars for spins as low as S = 0.03. 



has a drastically nonzero Lyapunov exponent. This ef- 
fect is referred to as a "chaos mimic" [26( and Hartl finds 
the same effect for knife edge orbits in the Kerr space- 
time @. This effect casts some doubt on the Lyapunov 
exponents calculated for the nonzero case. Some confi- 
dence is restored by the trend in the exponents as the 
spin decreases, specifically that the values approach zero 
as spin goes to zero. In Figure [13] we see a close up look 
at small spin values for the knife edge orbit. We still have 
the chaos mimic when S = 0, but we also see a very clear 
trend in the exponents as the spin decreases. 



Radial Distance 

FIG. 14: This is a plot ol the effective potential corresponding 
to a knife edge orbit constrained to high curvature regions. 
This orbit has L = 3.47 and E = 0.9432. The line E 2 = 
0.8897 is shown by the dotted line and shows that the particle 
is constrained to remain near r = 6m. 
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FIG. 15: Lyapunov exponents for spin values in the physical 
range. These orbits have orientation /3 = ir and constants of 
the motion for the spinless case are L = 3.47 and E — 0.9432. 
Note the chaos mimic when S = and the magnitude of the 
exponent is comparable to the larger values of Figure [7] 



C. Chaotic Orbits With Physical Spin 

We now present a particular chaotic orbit at spin values 
that may be relevant astrophysically. In this case, the 
constants of the motion are E = 0.9432 and L — 3.47 
and the spin vector orientation is /3 = n. These initial 
conditions define a knife edge orbit which is constrained 
to regions of high curvature. The effective potential for 
this orbit is shown in FigureQJJ Notice that the potential 
well confines the particle with E 2 = 0.8897, denoted by 
the dotted line, to remain close to r = 6m. 

In Figure [T5] we see Lyapunov exponents correspond- 
ing to this orbit. Similar to Figure [TU1 spin values greater 
than those shown on the graph cause the particle to be 
captured by the black hole. What is more interesting 
however, is that these nonzero Lyapunov exponents cor- 
respond to physical spin values and that these values have 
comparable magnitude to the larger exponents of Figure 

m 

One concern with this data is that the exponents seem 
to converge to the value of the chaos mimic of S = 0. 
We can use the KAM tori corresponding to this orbit to 
resolve this concern. In Figures [16] and [17] we compare 
the Poincare sections of this orbit when S = and equal 
steps between S = 3.0 x 10" 5 and S = 3.069 x 10~ 5 
respectively. 

Notice that when S = the phase space trajectory is 
confined to the surface of the 2 torus intersected by the 
r — P r plane. We can see this detail in the upper right 
of Figure [T7] In the nonzero spin cases the same Figure 
shows the breaking up of the tori surfaces which is indica- 
tive of chaos. Based on the positive indication for chaos 
given by both the Lyapunov exponent and the Poincare 
sections we conclude that this orbit is indeed chaotic for 



some initial conditions that correspond to astrophysically 
realizable amounts of spin. 

As the spin decreases from S = 3.069 x 10~ 5 the break- 
ing of the torus becomes less and disappears by 2.8 x 10~ 5 . 
Thus, it would appear that the points of Figure [T5l with 
smaller spins than this amount are not truly chaotic. 
Since these first points seem to converge to the chaos 
mimic of the S — case it is likely that the positive ex- 
ponents corresponding to unbroken tori are inflated by 
the chaos mimic exactly as in the zero spin case. We 
notice from the graph that the final point has smaller 
error bars than the previous points. Because this orbit 
is shown to be chaotic by its Poincare section, the true 
chaotic behavior of the particle may push it above the 
exponent magnitude created by the chaos mimic. 

VI. CONCLUSION 

We find that chaotic orbits with apparently astrophys- 
ically relevant spins do exist in the Schwarszchild space- 
time. In particular, we have shown that for a spin values 
between 5 = 3.0x 10~ 5 and S = 3.069 x 10~ 5 both the 
Lyapunov exponent and the Poincare section indicate a 
chaotic orbit. Recall that Hartl's @ bounds for chaotic 
spin values were 10~ 4 to 10~ 6 . These orbits combine 
the dynamics of knife edge orbits with the high curva- 
ture close to the center of the spacetime. While this is 
a special class of orbit, decaying orbits may exhibit this 
behavior before they are captured by the black hole. In- 
deed, Levin [23| has shown that compact binaries pass 
through a chaotic region when inspiraling and it is an 
interesting question as to whether decaying orbits could 
pass through these astrophysical, chaotic orbits. 
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The results of our analysis also put the cutoff spin value 
for chaotic orbits much lower than previously thought 
even for more general orbits. Suzuki and Maeda [I[ give 
a cutoff value of about S = 0.63 when the orbital angular 
momentum of the particle is L = 4. We have provided 
strong numerical evidence that chaotic orbits exist for 
spin values less than S = 0.01 for orbits with L = 3.6. 

Both [4| and 6] claim that the spin-spin interaction 
in the Kerr metric results in an even greater number of 
chaotic orbits with small spin than in the Schwarzschild 
spacetime considered here. Indeed, Hartl found a lower 
bound for spins that yield chaotic orbits of about S = 0.1. 
As the knife edge orbits that we have investigated here 
also exist in modified form in the Kerr spacetime, we 
might expect that there are even more chaotic orbits with 
physically relevant spins in Kerr. 



FIG. 16: These are Poincare sections in the r — P r plane for 
5 = and five orbits with spin magnitudes spaced equally 
from S = 3.0 x 10" 5 to S = 3.069 x 10" 5 . These orbits 
have the orientation /3 — ir and constants of the motion for 
the spinless case are L = 3.47 and E = 0.9432. Notice that 
unlike Figures[T]and[3]the initial conditions of the orbits make 
them difficult to distinguish at this level. 
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FIG. 17: This is a closer look at Fig 

is marked on the axis which spans r = 5.67 — 1 x 10~ D to 
r — 5.67 + 3.5 x 10~ 6 . This finer line in the upper left is the 
part of the section for the spinless case. Notice the distinct 
breaking of the tori for the other five orbits. These orbits 
increase in spin from top to bottom of the upper right corner. 



A final important question, which is the subject 
of 0, H| , is the potential impact of these chaotic orbits on 
the gravitational wave emission from these systems. Cer- 
tainly, the current results are suggestive that there is a 
chaotic regime in spin which might be astrophysically ac- 
cessible to extreme mass-ratio binaries. If so, then gravi- 
tational waves emitted by binary systems with large mass 
ratios might contain this chaotic imprint. Thus, analyz- 
ing these systems for gravitational wave emission would 
be a natural next step in understanding chaos in these 
systems. 
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